1. Field of the Invention
The present invention relates to a method for forming more rapidly a stochastic numerical model of Gaussian or Gaussian-related type, representative of the spatial distribution of a physical quantity (such as the permeability for example) in a porous heterogeneous medium (such as a hydrocarbon reservoir for example) calibrated in relation to data referred to as uncertain static and dynamic data.
2. Description of the Prior Art
Optimization in a stochastic context determines realizations of a stochastic model which meet a set of data observed in the field, referred to as static or dynamic data depending on the nature thereof. Static data correspond to observations on the studied physical quantity proper. When it is certain, the static datum is an exact value. In the opposite case, it is defined by a probability law. Dynamic data are characteristic of the displacement of fluids in the medium: they are, for example, production data (pressures obtained from well tests, flow rates, etc.).
In reservoir engineering, the realizations to be identified correspond to representations, in the reservoir field, of the distribution of carrying properties such as the permeability, the porosity or the facies distribution, each facies corresponding to a family of carrying properties. These realizations form numerical reservoir models. The available static data are, for example, punctual permeability, porosity or facies observations, and a spatial variability model determined according to punctual measurements. The punctual data can be defined by probability laws rather than exact values. For example, at a given point, a porosity value can be characterized by a normal probability law of mean 0.20 and variance 0.03. The dynamic data are directly related to the fluid flows in an underground reservoir, that is pressures, breakthrough times, flow rates, etc. The latter are often non-linearly related to the physical properties to be modeled. A randomly drawn realization is generally not in accordance with the whole of the data collected.
Static Data Integration
Coherence in relation to the static data is integrated in the model from kriging techniques
Journel, A. G., and Huijbregts, C. J., “Mining Geostatistics”, Academic Press, San Diego, Calif., 1978.
The general approach generates a non-conditional realization and in correcting it so that it meets the punctual observations and the spatial structure. Within this context, the punctual observations are assumed to be exact values (for example, at a given point, the permeability is 150 mD). Let there be a Gaussian and non-conditional realization y of the stochastic model Y. The corrected realization is obtained as follows:yc(x)=ydK(x)+[y(x)−yK(x)]
yc is the corrected and therefore conditional Gaussian realization. ydK and yK are obtained by kriging from the punctual observations and from the values of y simulated at the observation points.
When the punctual observations do not correspond to exact values, but are defined by probability laws, it is possible to use either a Bayesian approach, which is extremely time-consuming in making calculations, or a kriging technique. In the latter case, on which focus occurs herein, a preliminary stage is necessary to transform the probability law into a punctual value. This transformation cannot be any transformation: the punctual values obtained must meet the spatial structure and the probability laws from which they result. Clearly, the relation between probability laws and punctual values is not a single relation. Iterative transformation methods have been proposed for facies realizations, where the probability law is a uniform law, by:
Freulon, X., and Fouquet, C. de, “Conditioning a Gaussian Model with Inequalities”, in Geostatistics Troia '92, A. Soares, ed., Kluwer Academic Publishers, Dordrecht, The Netherlands, 1993,
Le Ravalec-Dupin, M., “Conditioning Truncated Gaussian Realizations to Static Data”, 8th Ann. Conf. Int. As. Math. Geol., Berlin, Germany, 15-20 Sep., 2002.
These techniques tend towards the desired probability density function, but very slowly. In practice, it is not desired to completely sample this function. A single set of punctual values meeting the spatial distribution and the probability laws is merely determined and it is used to obtain the corrected realization yc. It is this single and therefore invariable set that the optimization process constantly refers to so as to form realizations constrained simultaneously by the static and the dynamic data. This point is fundamental in the optimization process.
Dynamic Data Integration
Coherence in relation to the dynamic data is integrated in the model by means of an inverse procedure:
Tarantola, A., “Inverse Problem Theory—Methods for Data Fitting and Model Parameter Estimation”, Elsevier Science Publishers, 1987.
Recently, a geostatistical parameterization technique which simplifies the inverse problem has been introduced to constrain, by gradual deformation, the stochastic realizations to data on which they depend non-linearly. It forms the object of French Patents 2,780,798 and 2,795,841 filed by the assignee, and of the following publications, notably:
Hu, L. Y., 2000, Gradual Deformation and Iterative Calibration of Gaussian-Related Stochastic Models: Math. Geology, Vol. 32, No. 1,
Le Ravalec, M. et al., 2000, The FFT Moving Average (FFT-MA) Generator: An Efficient Numerical Method for Generating and Conditioning Gaussian Simulations: Math. Geology, Vol. 32, No. 6,
Hu, L. Y., Blanc, G. And Noetinger, B. (2001): Gradual Deformation and Iterative Calibration of Sequential Stochastic Simulations. Math. Geology, Vol. 33, No. 4.
This method has been successfully applied in various cases, notably from data obtained from oil development fields, as described in the following documents:
Roggero, F. et al., 1998, Gradual Deformation of continuous Geostatistical Models for History Matching, paper SPE 49004: Proc. SPE Annual Technical Conference and Exhibition, New Orleans,
Hu, L. Y. et al., 1998, Constraining a Reservoir Facies Model to Dynamic Data Using a Gradual Deformation Method, paper B-01: Proc. 6th European Conference on Mathematics of Oil Recovery (ECMOR VI), 8-11 Sep. 1998, Peebles, Scotland.
Assuming that fobs=(f1obs, f2obs, . . . , fMobs) represents all of the dynamic data collected in the field and f=(f1, f2, . . . , fM) the corresponding responses simulated for a realization yc already constrained by the static data are as explained above. In general, the responses f=(f1, f2, . . . , fM) are obtained by solving numerically the direct problem. Thus, if y0 represents a permeability field, data f can be pressure measurements. In this case, they are simulated from a flow simulator. The goal of a stochastic optimization is to produce realizations of Y which reduce the differences between the observed data and the numerically simulated corresponding responses. These differences are measured by the following objective function:
  J  =            1      2        ⁢                  ∑                  m          =          1                M            ⁢                                    ω            m                    ⁡                      (                                          f                m                            -                              f                m                obs                                      )                          2            
Coefficients ωm are weights assigned to data fmobs. fm are functions of realization y0 discretized over a large number of grid cells. In this sense, minimization of the objective function is a problem with several variables.
Let N be the number of grid cells forming realization y0. N is often very large (104˜107). It is therefore very difficult to perform an optimization directly in relation to the components of y0. Furthermore, realization y0, even modified, must remain a realization of Y. Parameterization by gradual deformation allows these difficulties to be overcome.
The gradual deformation technique allows construction of a continuous chain of realizations by combining an initial realization y0 of Y with another realization u1, referred to as complementary, of Y, u1 being independent of y0. The combination coefficients are for example cos(t) and sin(t), and the combined realization meets the relationy(t)=y0 cos t+u1 sin t where t is the deformation parameter.
Once the chain is formed, it can be explored by varying deformation parameter t and an attempt is made to identify, from among all the realizations of this chain, the realization which minimizes the objective function after integration of the static data by kriging. This minimization is performed in relation to t. Parameterization according to the gradual deformation method allows reduction of the number of dimensions of the problem from N to 1, where N is the number of values that constitute the field to be constrained. Furthermore, the sum of the combination coefficients squared being 1, the optimized realization still is a realization of Y: it follows the same spatial variability model as all the realizations of Y.
However, if the exploration of the realizations space is restricted to a single chain, the possibilities of sufficiently reducing the objective function are greatly limited. The above procedure therefore has to be repeated, but with new realization chains. These realization chains are constructed successively by combining an initial realization, which is here the optimum realization determined at the previous iteration, with a complementary realization of Y, randomly drawn each time. Thus, at iteration l, the continuous realization chain is written as follows:yl(t)=yl−1 cos t+ul sin t. 
yl−1 is the optimum realization defined at iteration l−1 and the ul are independent realizations of Y.
Minimizing the objective function in relation to t allows improvement or at least preservation of calibration of the data each time a new realization chain is explored. This iterative minimum search procedure is continued as long as data calibration is not satisfactory.
To date, the punctual values integrated by kriging in the model to account for the static data remain constant throughout the optimization process, even when they correspond to uncertain values. In the latter case, such a hypothesis can significantly slow down the optimization process and prevent minimization of the objective function.